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Abstract 

We apply a data assimilation techniques, inspired from meteorological 
applications, to perform an optimal reconstruction of the neutronic activ- 
ity field in a nuclear core. Both measurements, and information coming 
from a numerical model, are used. We first study the robustness of the 
method when the amount of measured information decreases. We then 
study the influence of the nature of the instruments and their spatial 
repartition on the efficiency of the field reconstruction. 

Keywords: Data assimilation, Schur complement, neutronic, activi- 
ties reconstruction, nuclear in-core measurements 



1 Introduction 

In this paper, we focus on the efficiency of a neutronic field reconstruction 
procedure with Data Assimilation when varying the number and the repartition 
of the available instruments. The data assimilation technique used for this 
reconstruction allows to combine, in an optimal and consistent way, information 
coming either from measurements or from a numerical model. 

Data assimilation methods are not commonly used in nuclear core physics [T] , 
contrary to meteorology or oceanography 2, 3, 4 . The procedure proposed here 
is the same as the one meteorologists use to obtain high accuracy meteorological 
reconstructed fields in time and space. This is the case, for example, of the 
commonly used meteorological re-analysis data set ERA-40 [S] among others 

mm- 

One of the main advantages of data assimilation is that it takes into account 
every kind of heterogeneous information within the same framework. Moreover, 
this method has a formalism that allows to adapt itself to instrument config- 
uration change. We exploit this last property here to study the quality of the 
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reconstructed activity field as a function of the number of available measure- 
ments. A major point in this study is to estimate the instrumented system 
robustness in the framework of a data assimilation reconstruction procedure. 
Moreover, such a study also informs about the effect of instrumentation design 
within a nuclear core and the resilience to instrument removal. 

In this paper, we first detail the data assimilation method and how it ad- 
dresses field reconstruction. To evaluate the influence of the number of in- 
struments on the activity field reconstruction, the repeated application of the 
method faces some huge computational issues. Those difficulties are overcomed 
using a matrix inversion method based on the Schur complement. A detailed 
presentation of this method is presented in Appendix|X] First we present the re- 
sults on a standard case with synthetic measurements and comments on them. 
To get a better understanding, we extend the results to other instrumental 
repartitions and other errors settings. This allows us to give some conclusions 
on the error and instrument repartition effects in activity field reconstruction 
using data assimilation. 

2 Data assimilation 

We briefly introduce the useful data assimilation key points to understand their 
use as applied in [5J [ID]- Data assimilation is a wider domain and these 
techniques are, for example, the keys of the nowadays meteorological operational 
forecasts [TT] . This is through advanced data assimilation methods that weather 
forecasting has been drastically improved during the last 30 years. All the 
available data, such as satellite measurements as well as sophisticated numerical 
models, are used. 

The ultimate goal of data assimilation methods is to estimate the inaccessible 
true value of the system state, x* where the t index stands for "true state" in 
the so called "control space". The basic idea is to combine information from an 
a priori on the state of the system (usually called x b , with b for "background"), 
and measurements (referenced as y°). The background is usually the result of 
numerical simulations, but can also be derived from any a priori knowledge. 
The result of data assimilation is called the analysis, denoted by x°, and it is 
an estimation of the true state x* we want to approximate. 

The control and observation spaces are not necessary the same, and a bridge 
between them needs to be built. This is the observation operator H, that trans- 
forms values from the space of the background to the space of observations. For 
our data assimilation purpose we will use its linearisation H around the back- 
ground. The inverse operation going from observation increments to background 
increments is given by the transpose H T of H. 

Two other ingredients are necessary. The first one is the covariance matrix 
of observation errors, defined as R = E[(y° — H(x t )).(y° — H(x t )) T ] where E[.] 
is the mathematical expectation. It can be obtained from the known errors on 
unbiased measurements which means E[y° — H(x. f )] — 0. The second one is the 
covariance matrix of background errors, defined as B = i?[(x & — x*).(x fc — x*) T ]. 
It represents the error on the a priori state, assuming it to be unbiaised following 
the E[x b — x'] = no biais property. There are many ways to get this a priori 
state and background error matrices. However, those matrix are commonly 
the output of a model and an evaluation of accuracy, or the result of expert 
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knowledge. 

It can be proved, within this formalism, that the Best Unbiased Linear 
Estimator (BLUE) x a , under the linear and static assumptions, is given by the 
following equation: 

x a = x 6 + K(y°- /fx 6 ), (1) 
where K is the gain matrix: 

K = BH T (HBH T + R) _1 . (2) 

Moreover, we can get the analysis error covariance matrix A, characterising the 
analysis errors x a — x*. This matrix can be expressed from K as: 

A = (I - KH)B, (3) 

where I is the identity matrix. 

It is worth noting that solving Equation [T] is, if the probability distribution 
is Gaussian, equivalent to minimise the following function J(x), x a being the 
optimal solution: 

J(x) = (x - x 6 ) T B- x (x - x") + (y° - HxjV 1 (y° - Hx). (4) 
This minimisation is known in data assimilation as 3D-Var methodology [8] . 

3 Data assimilation implementation 

The framework of this study is the standard configuration of a 900 MWe nu- 
clear Pressurized Water Reactor (PWR900) . To perform data assimilation, both 
simulation code and data are needed. For the simulation code, the EDF experi- 
mental calculation code for nuclear core COCAGNE in a standard configuration 
is used. The description of the basic features of this model are done in Section 

EEJ 

To have a good understanding of the instrumentation effect, we want to study 
various kind of configurations, even some that do not exist operationally and so 
cannot be tested experimentally. For that purpose, synthetic data are used that 
allows to have an homogeneous approach all along the document. Synthetic data 
is generated from a model simulation, filtered through an instrument model and 
noised according to a predefined measurement error density function (usually of 
Gaussian type). 

In the present case, we study the activity field reconstruction. An horizontal 
slice of a PWR900 core is represented on the Figure [TJ There is a total of 157 
assemblies within this core. Among those assemblies, 50 are instrumented with 
Mobiles Fissions Chambers (MFC). Those assemblies are divided verticaly in 
29 vertical levels. Thus, the size of the control x is 4553 (157 x 29). The size of 
the observation vector y° is 1450 (50 x 29). 

3.1 Brief description of the nuclear core model 

The aim of a neutronic code like COCAGNE is to evaluate the neutronic activity 
field and all associated values within the nuclear core. This field depend on the 
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Figure 1: The positions of MFC instruments in the nuclear core are localised 
in assemblies in black within the horizontal slice of the core. The assemblies 
without instrument are marked in white and the reflector, out of the reactive 
core, is in gray. 



position in the core and on the neutron energy. To do such an evaluation, 
the population of neutrons are divided in several groups of energy. In the 
present case only two groups are taken into account giving the neutronic flux 
<f> = ($i,$2) (even if the present code have no limit for the group number). 
The material properties depend on the position in the core, as the neutronic 
flux $, identified by solving two-group diffusion equations described by: 

-div(Digrad$i) + (E a i + S r )$i = - (vi?2 n $>i + ^2^/2*2) ^ 
-div(D 2 grad$ 2 ) + £ a 2$2 - S r $i = 0, 

where k is the effective neutron multiplication factor, all the quantities and the 
derivatives (except k) depend on the position in the core, 1 and 2 are the group 
indexes, S r is the scattering cross section from group 1 to group 2, and for 
each group, <£> is the neutron flux, E a is the absorption cross section, D is the 
diffusion coefficient, vllj is the corrected fission cross section. 

The cross sections also depend implicitly on the concentration of boron, 
which is a substance added in the water used for the primary circuit to control 
the neutronic fission reaction, throught a feedback supplementary model. This 
model takes into account the temperature of the materials and of the neutron 
moderator, given by external thermal and thermo-hydraulic models. A detailed 
description of the core physic and numerical solving can be found in reference 
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The overall numerical resolution consists in searching for boron concentration 
such that the eigenvalue k is equal to 1, which means that the nuclear power 
production is stable and self-sustaining. It is named critical boron concentration 
computation. 

The activity in the core is obtained through a combination of the fluxes $ = 
($1, $2), given on the chosen mesh of the core. Using homogeneous materials for 
each assembly (for example 157 in a PWR900 reactor), and choosing a vertical 
mesh compatible with the core (usually 29 vertical levels), this result in a field 
of activity of size 157 x 29 = 4553 that cover all the core. 

3.2 The observation operator H 

The H observation operator is the composition of a selection and of a normali- 
sation procedure. The selection procedure extracts the values corresponding to 
effective measurement among the values of the model space. The normalisation 
procedure is a scaling of the value with respect to the geometry and power of 
the core. The overall operation is non linear. However, with a range of value 
compatible with assimilation procedure, we can calculate the linear associated 
operator H. This observation matrix is a (4553 x 1450) matrix. 

3.3 The background error covariance matrix B 

The B matrix represents the covariance between the spatial errors for the back- 
ground. In order to get those, we estimate them as the product of a correlation 
matrix C by a normalisation factor. 

The correlation C matrix is built using a positive function that defines the 
correlations between instruments with respect to a pseudo-distance in model 
space. Positive functions have the property (via Bochner theorem) to build 
symmetric defined positive matrix when they are used as matrix generator [131 
114"] . In the present case, Second Order Auto- Regressive (SOAR) function is 
used to prescribe the C matrix. In such a function, the amount of correlation 
depends from the euclidean distance between spatial points. The radial and 
vertical correlation length (L r and L z respectively, associated to the radial r 
coordinate and the vertical z coordinate) have different values, which means we 
are dealing with a global pseudo euclidean distance. The used function can be 
expressed as follow: 

The matrix obtained by the above Equation [6] is a correlation matrix. It is 
then multiplied by a suitable variance coefficient to get covariance matrix. This 
coefficient is obtained by statistical study of difference between model and mea- 
surements in real case. In our case, the size of the B matrix is related to the 
size of model space so it is (4553 x 4553). 

3.4 The observation error covariance matrix R 

The observation error covariance matrix R is approximated by a diagonal ma- 
trix. This means it is assumed that no significant correlation exists between the 
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measurement errors of the MFC. The usual modelling is to take those value as 
a percentage of the observation. This can be expressed as: 

R « = ( a (y°)]f . Vj (7) 

The parameter a is fixed according to the accuracy of the measurement and the 
representative error associated to the instrument. The size of the R matrix is 
related to the size of observation space, so it is (1450 x 1450). 



4 General results on instrument removal 

To test the robustness, many BLUE calculations need to be done to evaluate 
the results quality with instruments configuration modifications. We want to 
have an evaluation of the quality of reconstruction as a function of the num- 
ber of instruments, with a significant statistical result. To efficiently perform 
these numerous computations, a specific method using Schur complement was 
developped. The details of this new method are reported in Appendix [Xj 

Here, we are interested in the evaluation of the quality of the analysis x a as a 
function of the amount of provided information. To quantify this effect we make 
a statistic of 200 scenarios of instruments removal. We are making those statis- 
tics on several hypothesis, starting from a complete instrument configuration 
and then removing instruments two by two until none remains. The calcula- 
tion are done on the basis of the algorithm and hypothesis on data assimilation 
described previously. 

To quantify the impact of removed instruments on the analysis, we look at 
the percentage quantity v defined as follow: 

\\y° ef - Hx b \\ - \\y° ef -Hx a \\ 

where x° e j corresponds to the analysis when no instrument is removed (this is 
the best estimation possible with respect to the information available on the 
system), and where y° e f and H are the reference observations and observation 
operator used to build x° e j. H stand for the observation operator when no 
instrument are removed. This criterium, which is basd on the norm of the 
innovation vector y° e f — Hx b , focuses on measurements. Since ||y° e y — -Hx a || 
is greater than \ \y° e f — Hx^ e f\\ (best estimate) and smaller than ||y° e/ — -ffx b || 
(innovation), v is a measure of the quality of the analysis. 

Such a definition have several advantages. First of all, the limit of this 
function are interesting. On the one hand, limit when no instrument is removed 
is 100%. On the other hand, the limit when all instruments are removed is 
0%. With such a formula we can compare the variation of the information on 
a unique scale. If we obtain some value above the limit of 100%, this mean the 
parametrisation of data assimilation was not done correctly. 

The interest of using this formula is that it can be applied directly as well 
to experimental data as to synthetic data without any change. 

On Figure [2] are presented the results of the quantity v as a function of the 
number of removed instruments. 

As expected, the relative quality of reconstruction decreases as a function of 
the number of instruments removed. However within this decrease, three phases 
can be seen: 
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Efficiency of reconstruction as a function of instument removal 
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linear fit slope -1 .64 
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Figure 2: The plain curve represents the value v given by Equation [8] as a func- 
tion of the number of instruments removed. This results come from a typical 
PWR900 instrument location. The two other lines are linear regression corre- 
sponding to both decreasing regime. 



1. A first phase of slow decreasing until we removed roughly 20 instruments. 
This phase is rather clear and can be fitted by a linear regression with a 
slope of —1.64 (arbitrary unit (a.u.) per instrument). The fit is shown in 
(green) dash line on Figure [2j 

2. After 20 instrument removed the decreasing speed of the slope increases. 
The second linear fit has a slope of —2.28 (a.u. per instrument). This fit 
is shown in (blue) dotted line in Figure [2j 

3. Beyond 40 instruments removed, we reach a third phase of stagnation then 
a brutal decrease to the limit value imposed by Equation [8j 

This characteristic behaviour can be seen in several cases that we studied. 
We have also noticed it on real measurements [T5J [TB] . The transition between 
the two first decreasing phases is specially strong when we do the analysis using 
real measurements. 

First, we explain why the third phase is marked by a stagnation of the mean 
value of ||y° — Hx a || over the set of removal scenarios taken into account. To 
understand that effect, we work on both the cases where two instruments are 
removed (i.e 48 are remain), and the ones where two instruments are remaining. 
On those two cases, at most 1225 scenarios are possibles, corresponding to the 
C|° combinations. Using the hybrid Schur method, we can calculate all those 
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cases with a rather cheap computing time to obtain a good statistics. We plotted 
the distribution of the value of ||y° — Hx a || over all the scenarios in Figure [3j 
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Figure 3: Distribution of the norms of ||y° — Hx a || for all the cases where only 
two instruments are suppressed or two instruments remain in the instrument 
network on a PWR900 reactor. 



On the Figure [3] we notice a big difference between the two distributions 
plotted, that correspond to the cases where two instruments are removed (i.e. 
48 remain) or two remain. 

The shifting of the mean value between those two extreme cases is logical as 
available information is dramatically changing. However, the shape of the dis- 
tribution is also vastly changing. We move from a very sharp distribution, when 
48 instruments remain, to a rather broad one when two instruments remain. 

The first interesting point of this shape change is that all the instruments do 
not have the same influence on the activity field reconstruction. To understand 
better this effect let's assume that all the instruments are equivalent. In this 
case, as the number of scenarios present in the two distributions of Figure [3] 
are the same, only a shifting of the mean value should have be seen. However 
we have not only a translation of the distribution but also a broadening. Thus, 
there is a non equivalence of instruments within the data assimilation procedure, 
in terms of marginal information assigned to each instrument, depending on its 
location in the core 

The second point of interest is that the broadening of the distribution is 
asymmetric. The distribution extends towards the higher values of norm. This 
effect explains the stagnation of the v quantity when few instruments remain. 
The source of stagnation is the discrepancy between the most probable value 
of the distribution and the mean value of the distribution which is higher. The 
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most probable value leads to a decrease without stagnation. However, looking at 
the mean value (that has more physical meaning), we see that this one stagnates 
due to the asymmetric broadening that compensates the overall decrease in 
mean value. 

Now the origin of the two slopes in decreasing of the information represented 
by the Figure [2] will be investigated. The repartition of the instruments in a 
standard PWR900 is very complex as shown in Figure [l] This complexity of 
the repartition does not make the situation easy to understand. Thus, we want 
to study the case on a simpler repartition of the instrumentation to see if the 
effect of two phases decreasing persists. 

5 Repartition effects in instrument removal 

The position of the instruments are presented in Figure |4j The core geome- 
try and assemblies configuration is the same as a PWR900 one, however the 
instrument are located regularly on a Cartesian map. 




2 4 6 8 10 12 14 

x position 



Figure 4: The MFC instruments within the nuclear core are localised in assem- 
blies in black within the horizontal slice of the core. The assemblies without 
instrument are marked in white and the reflector is in gray. 

Within this configuration, only 40 MFC are used, which is a bit less than 50 
of the standard PWR configuration presented on Figure[TJ With this repartition, 
we do the same analysis as on the previous one. The evolution on v as a function 
of the number of removed instruments is plotted in Figure [5j 

In the Figure [5j the quality of the activity reconstruction decreases quasi- 
linearly as a function of the number of removed instruments. This goes on until 
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Efficiency of reconstruction as a function of instument loss 
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Figure 5: The plain curve represents the value v given by Equation [8] as a 
function of the number of instruments removed. This results come from a regular 
instruments location in a PWR900. The dashed line represents the linear fit of 
the steady part of the curve. 



we reach the stagnation phase. It appears, comparing Figure [5] to [2] that the 
variation of the decrease slope is related to the geometrical repartition of the 
instruments. Such uniform linear decreasing can be observed also in several 
other rather geometrically regular repartition, as one with repartition on diago- 
nal line. Some of the repartitions have densities of instruments close to the one 
of the standard PWR900, which do not change the overall behaviour. 

Looking at the slope factor of the linear fit, we notice that the one with 
regular MFC repartition (Figure [5 ) is in the range of the slopes obtained with 
standard MFC repartition (Figure 2l). 

Thus the removal or fault resilience of the instrument set seems to depend 
on the transition point between the two decreasing steps. In the lower range of 
instrumental density, it is better to have a regular repartition, and in the higher 
one, it is better to have a complex ad-hoc repartition. In real PWR900 nuclear 
core, because the set of intruments is fixed at a high density level, these results 
indicate that it is more robust to have an ad-hoc repartition of the instruments 
as for now. 



6 Conclusion 

We proposed and studied here an original method to test how the neutronic 
activity field reconstruction by data assimilation is tolerant to information re- 
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moval. The core of the reconstruction method is based on data assimilation, 
that is widely use in earth sciences, and allows a very good reconstruction of the 
activity all over the nuclear core. An hybrid method for fast matrix inversion 
partially based on Schur complement allows the execution on numerous analy- 
sis from which statistical results are derived. For all these analyses, synthetic 
data are used to try non-experimental intruments repartitions, but similar re- 
sults were established using real data. This application on real data prove both 
the reliability and the quality of the calculation code and the data assimilation 
methodology. 

Using those advanced calculation methods, it was shown that the slopes of 
the reconstruction quality is mainly governed by repartition for the instruments. 
Depending on the chosen repartition, the decrease consists in two or three dis- 
tinct phases. The ultimate stagnation phase in this decreasing is governed by 
both statistical effect and heterogeneity of instruments influence. 

The behaviour with two phases within the decreasing quality of the recon- 
struction as a function of the number of instruments removed is understood in 
term of repartition effect, but not quantified. However it can be seen as a phase 
transition between two states of instrumental configuration. The quantification 
of this transition is worth studying. 
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A Appendix: Schur complement method to op- 
timise calculation 

Within the BLUE assimilation method, the limiting factor in calculation time 
is the matrix inversions. In Equation [2j the costly part is the inversion of the 
term: 

M = HBH T + R. (9) 

The inversion cost on huge matrix as M (around 4000 x 4000 in the present 
case) was such that the time calculation of the above evaluation was extremely 
time consuming. Then we had to optimise the computing cost. 

We noticed that the calculations are more time consuming when only few 
instruments are removed. In this case the M matrix is still huge. 

Thus, the idea is to use the information obtained in the inversion of the full 
size matrix to shorten calculation, to calculate smaller size matrix in a reason- 
able time. In this case, we want to calculate the new matrix as a perturbation 
of the original one. Such a method exists and exploits the Schur complement of 
the matrix. 

We assume we want to suppress some instruments to a given configuration. 
With respect to the Equation [2j we need to calculate a new matrix K„. The n 
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index is standing for referring the new matrix we want to calculate. For that 
according to Equation [9] we have to determine a new matrix M„. 

This determination of M„ is obtained from the knowledge of the invert of 
the matrix M g calculated over all the instruments. The indices g is used to 
denote the reference matrix we start from M g according to Equation [9] 

All the components of the new matrix K„ can be obtained by suppressing 
the lines and columns corresponding to removed instruments in M fl , inverting it 
and then multiplying this matrix by the corresponding H„ and B„. We notice 
that, in our case, we get B„ = B 9 as we do not affect the model space. 

To make the demonstration easier, but without losing any generality, we can 
assume that the suppressed instruments correspond to the lower square of M 9 . 
If it is not the case, it is always possible to reorganise matrix in such a way. 

Now we put the M g under a convenient form, separating remaining measures 
from removed ones. Assuming the starting matrix M 9 is m x m and that we 
plan so suppress s measurements, we can write M g in the following way: 

where: 

• P g contains the remaining measurements, and is a p x p matrix, 

• S g contains the suppressed measurement, and is a s x s matrix, 

• Q g ct R g represent the dependence between remaining measured and sup- 
pressed ones. In the particular case we are dealing with, we have Q J = R s . 
However, no further use of this property is done. 

With such a decomposition, we got the equality m = p + s. 

The P g matrix corresponds to the remaining instruments, thus we have the 
equality: 

P S = M„, (11) 

The decomposition given in Equation [lO] is the one required to build the Schur 
complement of this matrix [17]. Under the condition that P g can be inverted, 
the Schur complement is the following quantity: 

S g - RgP^Qg, (12) 

and is noted (M 9 /P 5 ). This notation reads as Schur complement of M 5 by P g . 

Thus we look for a cheap way to calculate P~ 1 knowing M" 1 . For that, we 
use the Banachiewiz formula [17] that gives invert of M 5 as a function of P 9 , 
Q g , Tig, Sg SLiid (M g /P g ) matrices: 

M ^ - ill sfV 1 (13) 



"9 ^9 

Pg" 1 + Pg- 1 Qg(Mg/Pg)^ 1 RgPgT 1 - P g ^ Q g ( Mg / P g , 

yp g )-^R g p-- 



-(Mg/Pg)- 1 RgP- 1 (MjP^ 
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We define the 4 sub-matrices P g , Q g , R g and S g by: 



P g = P S - 1 +P J - 1 Q ( (M S /P S )- 1 H ! P; 1 > (14) 

Q g = -P^QgiMg/Pgy 1 , (l5) 

R g = -(rvVPj-^P; 1 , (16) 

Sg - (Mg/Pg)- 1 . (17) 

Rearranging those terms we get: 

p; 1 = p 9 - <Wr 5 . (18) 



As by hypothesis we know the inverse M g 1 of the global matrix, we are able to 
extract P g , Q g , R g and S g from the whole inverted matrix. Thus, the main cost 
to obtain the inverse of P g of size p x p becomes the one of inverting S g which 
size is q x q. In first approximation, if the number of measurements to suppress 
is smaller than the number of remaining one, this methods gives a notable gain. 
As soon as the matrix P g 1 = M^ 1 , final calculation of K„ is straightforward. 

To highlight the advantages of this method with respect to the standard 
inversion of sub-matrix, some tests are shown on a 4000 x 4000 full semi-definite 
positive regular matrix. The curves showing the effective computing time in 
percentage of the computing time of the full matrix are presented in Figure [6| 

Figure [6] shows that, when sub-matrix has roughly the size of the initial 
matrix, the inversion by Schur complement is far more efficient than the direct 
inversion. Above around 60% of the size of the initial matrix, the direct inver- 
sion becomes more efficient. The crossing point is at 60% instead of 50%, as 
expected in first approximation. This difference comes from the few additional 
multiplications that need to be done in the Schur complement calculation, as 
we can see in Equation [12] Globally, we see that the most efficient method is 
to use an hybrid calculation that choose the best way to make the calculation 
as a function of the number of measurements removed. To quantify the im- 
provement of such an hybrid choice, we integrate the curves of direct inversion 
and compare it to the integral of hybrid option (minimum of both cases) within 
the instrument loss range. The ratio in percentage of both integrals shows that 
benefit of this hybrid method represents an overall gain of 64% with respect to 
standard method. 
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Figure 6: Calculation time as a function of the size of the sub-matrix. Calcula- 
tion time is given as a percentage of the calculation time of the full size matrix. 
The size of the sub-matrix is given a a percentage of the size of the full matrix. 
The curves in full line (red) and in dashed line (green) represent the inversion 
by Schur complement and direct inversion respectively. 
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